#include "pgroup.h"
#include "main.h"


pgroup::pgroup(vector<vec3> R)
{
	r = R;
	rcell = 3.5;
	initialize0();
}




pgroup::pgroup(vector<vec3> R, real rcell0)
{
	r = R;
	rcell = rcell0;
	initialize0();
}




void pgroup::initialize0()
{
	isx    = int(XSIZE/rcell);
	isy    = int(YSIZE/rcell);
	isz    = int(ZSIZE/rcell);
	rxsize = real(isx)/XSIZE;
	rysize = real(isy)/YSIZE;
	rzsize = real(isz)/ZSIZE;
	if(isx<3 || isy<3 || isz<3)
	{
		cout<<"system size is too small";
		exit(-1);
	}	
	iside  = isx;
	iside2 = isx * isy;
	iside4 = iside + iside2;
	nbox   = iside2*isz;
	ibox.resize(nbox);
	_debug = false;
}




void pgroup::initialize()
{
//	npair = p.size();
	ccc.assign(npair,1.);
	ddd.assign(npair,1.);
	rc.assign(npair,2.5);
	rlist.assign(npair,3.5);
	vlist.clear();
	vlist.resize(npair);
}




int pgroup::pair_id(int d1, int d2)
{
	return 0;
}




void pgroup::labels()
{	
	int nc0, nc1, nc2, nc3, nc4, nc5, nc6, nc7, nc8, nc9, nc10, nc11, nc12, nc13;
	int ip, jp, kp, jr, kr;
    unsigned int i, j, k;
		
	for (i = 1; i <= isx; i++) {
		for (j = 1; j <= isy; j++) {
			for (k = 1; i <= isz; k++) {
				ip   = (i == isx ? 1 : i+1);
				jp   = (j == isy ? 1 : j+1);
				kp   = (k == isz ? 1 : k+1);
				jr   = (j == 1 ? isy : j-1);
				kr   = (k == 1 ? isz : k-1);
				
				nc0  = i  +  j*iside +  k*iside2 - iside4;
				
				nc1  = ip +  j*iside +  k*iside2 - iside4;
				nc2  = i  + jp*iside +  k*iside2 - iside4;
				nc3  = i  +  j*iside + kp*iside2 - iside4;
				nc4  = ip + jp*iside +  k*iside2 - iside4;
				nc5  = ip +  j*iside + kp*iside2 - iside4;
				nc6  = i  + jp*iside + kp*iside2 - iside4;
				nc7  = ip + jp*iside + kp*iside2 - iside4;
				nc8  = ip + jr*iside +  k*iside2 - iside4;
				nc9  = ip +  j*iside + kr*iside2 - iside4;
				nc10 = ip + jr*iside + kr*iside2 - iside4;
				nc11 = i  + jp*iside + kr*iside2 - iside4;
				nc12 = ip + jr*iside + kp*iside2 - iside4;
				nc13 = ip + jp*iside + kr*iside2 - iside4;
				
				label.push_back(make_pair(nc0,nc1));
				label.push_back(make_pair(nc0,nc2));
				label.push_back(make_pair(nc0,nc3));
				label.push_back(make_pair(nc0,nc4));
				label.push_back(make_pair(nc0,nc5));
				label.push_back(make_pair(nc0,nc6));
				label.push_back(make_pair(nc0,nc7));
				label.push_back(make_pair(nc0,nc8));
				label.push_back(make_pair(nc0,nc9));
				label.push_back(make_pair(nc0,nc10));
				label.push_back(make_pair(nc0,nc11));
				label.push_back(make_pair(nc0,nc12));
				label.push_back(make_pair(nc0,nc13));
				
				// periodic boundary condition need to be considered here
				// more codes will be added here

			}
		}
	}
	
	nlabel = label.size();
}




void pgroup::verlet()
{
	atom temp_atom;
	vec3 vatom;
	int abox, value;
    unsigned int i, i1, i2, id1, id2;
	real dr;
	
	ibox.clear();
	ibox.resize(nlabel);				// clear ibox
	
	for (i = 1; i<r.size(); i++) {
		abox = int(r[i].x*rxsize+1.)+int(r[i].y*rysize+1.)*iside+int(r[i].z*rzsize+1.)*iside2-iside4;
		ibox[abox].push_back(i);
	}									// renew ibox
	
	vlist.clear();
	vlist.resize(npair);				// clear vlist
	
	for (i = 0; i<nbox; i++) {		// build verlet list from same cell
		vector<int>& t_list=ibox[i];
		
		if(t_list.empty())continue;
		
		for (i1 = 0; i1<t_list.size()-1; i1++) {
			id1 = t_list[i1];
			for (i2 = i1+1; i2<t_list.size(); i2++) {
				id2   = t_list[i2];
				value = pair_id(id1,id2);
				if(value){
					dr = (r[id1] - r[id2]).module();
					if(dr<rlist[value]){
						vlist[value].push_back(make_pair(id1,id2));
					}
				}
			}
		}
	}
	
	for (i = 0; i<nlabel; i++) {		// build verlet list from neigbor cells
		vector<int>& t_list1=ibox[label[i].first];
		vector<int>& t_list2=ibox[label[i].second];
		
		if(t_list1.empty() || t_list2.empty())continue;
		
		for (i1 = 0; i1<t_list1.size(); i1++) {
			id1 = t_list1[i1];
			for (i2 = 0; i2<t_list2.size(); i2++) {
				id2   = t_list2[i2];
				value = pair_id(id1,id2);
				if(value){
					dr = (r[id1] - r[id2]).module();
					if(dr<rlist[value]){
						vlist[value].push_back(make_pair(id1,id2));
					}
				}				
			}
		}
	}
}




void pgroup::eval()
{
	
}




void pgroup::debug()
{
	if(_debug){
		cout<<"warning: debug mode already turned on"<<endl;
	}
	else {
		_debug = true;
		cout<<"**********Debug Mode**********"<<endl;
	}
	// more codes shall be added here
}
